A map chart is used to show items on a background that is often geographical.
Map charts in R can contain different type of information in one or more different layers. Maps can contain interactive shapes or display markers of different types on an image or map background.
Steps to draw a map using R: 1). Find the required map data, such as latitude and longitude, boundary information, etc., and load it into R; 2). Data reconstruction meets the requirements of function to draw the map.
The package maps comes with maps of the world, maps of Italy, Canada, the United States, France, New Zealand and other regions. But the available maps are limited by the R package.
library(maps)
library(ggmap)
## Loading required package: ggplot2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(leaflet)
library(htmlwidgets) # replacing leaflets
library(RColorBrewer)
The ggmap package can link to Google Map, OpenStreetMap and Stamen Maps, and provide latitude and longitude, address, distance and other related information. The package ggmap needs to register to get the API key first.
map("world", fill = TRUE, col = 'lightBlue',
ylim = c(-60, 90), mar = c(0, 0, 0, 0))
library(RColorBrewer)
map('state', region = c('california', 'oregon', 'washington'),
fill = TRUE, col = colorRampPalette(brewer.pal(8,'Pastel1'))(3), mar = c(2, 3, 4, 3))
You need to specify the coordinates to draw a certain area and zoom in
leaflet() %>%
addTiles() %>%
setView( lng = -95.09,
lat = 29.55,
zoom = 10) %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") ## operated by the NASA data
A choropleth map is a type of thematic map in which areas are shaded or patterned in proportion to a statistical variable that represents an aggregate summary of a geographic characteristic within each area, such as population density or per-capita income.
Choropleth maps provide an easy way to visualize how a measurement varies across a geographic area or show the level of variability within a region. A heat map or isarithmic map is similar but does not use a priori geographic areas. They are the most common type of thematic map because published statistical data (from government or other sources) is generally aggregated into well-known geographic units, such as countries, states, provinces, and counties.
library(RColorBrewer)
library(rgdal)
## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.3, released 2021/04/27
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgdal/proj
library(ggplot2)
library(ggmap)
library(viridis)
## Loading required package: viridisLite
##
## Attaching package: 'viridis'
## The following object is masked from 'package:maps':
##
## unemp
library(leaflet)
library(htmltools)
library(mapproj)
library(usmap)
library(ggplot2)
library(data.table)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1
## ✓ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::between() masks data.table::between()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks data.table::first()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks data.table::last()
## x purrr::map() masks maps::map()
## x purrr::transpose() masks data.table::transpose()
library(geofacet)
head(statepop)
## # A tibble: 6 × 4
## fips abbr full pop_2015
## <chr> <chr> <chr> <dbl>
## 1 01 AL Alabama 4858979
## 2 02 AK Alaska 738432
## 3 04 AZ Arizona 6828065
## 4 05 AR Arkansas 2978204
## 5 06 CA California 39144818
## 6 08 CO Colorado 5456574
plot_usmap(data = statepop, values = "pop_2015", color = "white") +
scale_fill_continuous(name = "Population 2015", label = scales::comma) +
theme(legend.position = "right")
plot_usmap(data = countypop, values = "pop_2015", color = "white") +
scale_fill_continuous(name = "Population (2015)", label = scales::comma) +
theme(legend.position = "right")
We can even use map and ggplot together to built up the plot
library(tidyverse)
library(ggplot2)
library(maps)
library(usmap)
state = map_data("state")
head(state)
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
## 4 -87.53076 30.33239 1 4 alabama <NA>
## 5 -87.57087 30.32665 1 5 alabama <NA>
## 6 -87.58806 30.32665 1 6 alabama <NA>
state$region=tolower(state$region)
statepop$full=tolower(statepop$full)
state_pop <- merge(as_tibble(state), statepop, by.x = "region", by.y="full")
p <- ggplot() +
geom_polygon( data=state_pop,
aes(x=long, y=lat, group=group,
fill = pop_2015/1000000),
color="white", size = 0.2) +
scale_fill_continuous(name="State Population",
low = "#ffe6cc",
high = "#ffa700",
limits = c(0,40),
breaks=c(5,10,15,20,25,30,35),
na.value = "grey50")
p
The cartogram heatmap does not express the value by changing the shape or size of the map like cartogram. Each area cartogram heatmap in is represented by a square of the same size, and the value is expressed by color. It is widely used in the data of the states in the United States than other regions.
if(!require(statebins)) install.packages("statebins")
## Loading required package: statebins
library(statebins)
library(tidyverse)
library(viridis)
covid <- read.csv("covid_states.csv")
head(covid)
## submission_date state tot_cases conf_cases prob_cases new_case pnew_case
## 1 02/12/2021 UT 359641 359641 0 1060 0
## 2 03/01/2021 CO 438745 411869 26876 677 60
## 3 08/22/2020 AR 56199 NA NA 547 0
## 4 07/17/2020 MP 37 37 0 1 0
## 5 08/12/2020 AS 0 NA NA 0 0
## 6 05/22/2021 MA 704796 659246 45550 451 46
## tot_death conf_death prob_death new_death pnew_death created_at
## 1 1785 1729 56 11 2 02/13/2021 02:50:08 PM
## 2 5952 5218 734 1 0 03/01/2021 12:00:00 AM
## 3 674 NA NA 11 0 08/23/2020 02:15:28 PM
## 4 2 2 0 0 0 07/19/2020 12:00:00 AM
## 5 0 NA NA 0 0 08/13/2020 02:12:28 PM
## 6 17818 17458 360 5 0 05/23/2021 01:37:59 PM
## consent_cases consent_deaths
## 1 Agree Agree
## 2 Agree Agree
## 3 Not agree Not agree
## 4 Agree Agree
## 5
## 6 Agree Agree
#preparation
covid$state<-as.character(covid$state)
covid$tot_death<-as.numeric(covid$tot_death)
#draw the plot
statebins(covid, #dataframe
state_col="state", #the name of state
value_col = "tot_death", #the number of death
round=TRUE,
dark_label = "white", #dark label as white
light_label = "black")+ #light label as black
scale_fill_viridis(option = "E",begin = 0.95,end = 0,name="Deaths")+ #the color of fill
labs(title = "COVID-19 U.S. Cumulative Deaths by December 12,2021", #title and caption
caption="Source: Centers for Disease Control and Prevention")+
theme_statebins(legend_position = "right")
We could seperate the total deaths into sub-groups
##preparation
covid$tot_death[which(covid$tot_death<=100)]<-1
#group the number of tot_death
covid$tot_death[which(covid$tot_death>100&covid$tot_death<=1000)]<-2
covid$tot_death[which(covid$tot_death>1000&covid$tot_death<=5000)]<-3
covid$tot_death[which(covid$tot_death>5000&covid$tot_death<=10000)]<-4
covid$tot_death[which(covid$tot_death>10000&covid$tot_death<=30000)]<-5
covid$tot_death[which(covid$tot_death>30000)]<-6
covid$tot_death<-as.factor(covid$tot_death) # change them to factor
ggplot(covid,aes(state = state, fill=tot_death))+
geom_statebins()+
scale_fill_brewer(palette = 3,
direction = -1,
name="Cumalative Deaths",
limits=c("6","5","4","3","2","1"), #order the group
labels=c(">30000","10001-30000",
"5001-10000","1001-5000",
"101-1000","<=100")) + # chage the label
labs(title = "COVID-19 U.S. Cumulative Deaths by December 12, 2021", #add title
caption="Source: Centers for Disease Control and Prevention")+ #caption
theme_statebins(legend_position = "right")
A facet map is faceting the geographic regions and draws the data according to the corresponding geographic location. The graph can intuitively compare the distribution of data in different regions.
Package: Facet map uses the geofacet package to perform faceting according to different geographic areas, and make the relative position of each area consistent with the actual geographic location.
library(ggplot2)
library(geofacet)
Now I am using the election dataset in the geofacet which is the results from 2016 election
head(election)
## state candidate votes pct
## 1 Alabama Clinton 729547 34.357946
## 2 Alabama Trump 1318255 62.083092
## 3 Alabama Other 75570 3.558962
## 4 Alaska Clinton 116454 36.550871
## 5 Alaska Trump 163387 51.281512
## 6 Alaska Other 38767 12.167617
ggplot(election,aes(candidate,pct,fill=candidate))+
geom_bar(stat = 'identity')+
scale_y_continuous(limits = c(0,100),breaks = seq(0,100,30),minor_breaks = T)+
scale_fill_manual(values = c("#4f6f57", "#99ccb3", "#03c03c"))+
labs(title = '2016 US Election')+
facet_geo(~state,grid =us_state_grid1 )+ ## facet by the geographic regions
theme_bw()+
theme(
axis.text.x = element_blank()
)
## Note: You provided a user-specified grid. If this is a generally-useful
## grid, please consider submitting it to become a part of the geofacet
## package. You can do this easily by calling:
## grid_submit(__grid_df_name__)
A connection map shows the connections between several positions on a map. The link between 2 locations is usually drawn using great circle: the shortest route between them. We are going to use the geosphere package.
library(ggplot2)
library(maps)
library(geosphere)
##
## Attaching package: 'geosphere'
## The following object is masked from 'package:htmltools':
##
## span
airports <- read.csv("airports.csv",header = T)
flights <- read.csv("flights.csv",header = T,as.is=TRUE)
dataVI <- read.csv("thebeltandroad.csv",header = T)
dataVI
## loc long lat type
## 1 Beijng 116.397128 39.916527 1
## 2 Xian 108.934250 34.230530 1
## 3 Lanzhou 103.718780 36.103960 1
## 4 Urumqi 88.311040 43.363780 1
## 5 Kazakhstan 71.465669 51.161882 1
## 6 Kyrgyzstan 74.567320 43.103590 1
## 7 Tajik 68.787038 38.568123 1
## 8 Uzbekistan 69.248888 41.314031 1
## 9 Iran 51.384950 35.691309 1
## 10 Turkey 32.856293 39.936684 1
## 11 Russia 37.615010 55.757770 1
## 12 Germany 13.403229 52.523688 1
## 13 Netherlands 4.885376 52.375914 1
## 14 Italy 12.478860 41.953351 1
## 15 Greece 23.701211 37.991796 2
## 16 Kenya 36.800086 -1.286858 2
## 17 India 78.083620 21.175326 2
## 18 Sri Lanka 80.622447 7.654822 2
## 19 Singapore 103.784644 1.287574 2
## 20 Malaysia 102.092093 4.321244 2
## 21 Philippines 121.354063 12.900646 2
## 22 Vietnam 105.697963 21.002657 2
## 23 Beihai 109.119860 21.502835 2
## 24 Haikou 110.334082 19.979834 2
## 25 Guangzhou 113.273050 23.143924 2
## 26 Quanzhou 118.663458 24.897163 2
## 27 Fuzhou 119.298165 26.082742 2
Mark the location for airports
data1 <- dataVI[dataVI$type == "1",] #sea and land
data2 <- dataVI[dataVI$type == "2",]
maps::map('world',
col="grey", fill=TRUE, bg="white", lwd=0.05,
mar=rep(0,4),
border=0,
xlim=c(-20, 150), ylim=c(-20, 70)
)
points(x=data1$long, y=data1$lat, col="tomato", cex=1, pch=20) #cex: size of point,pch: shape of the point
points(x=data2$long, y=data2$lat, col="slateblue", cex=1, pch=20)
maps::map('world',
col="darkseagreen1", fill=TRUE, bg="white", lwd=0.05,
mar=rep(0,4),
border=0,
xlim=c(-20, 150), ylim=c(-20, 70)
)
points(x=data1$long, y=data1$lat, col="tomato", cex=1, pch=20) #cex: size of point,pch: shape of the point
points(x=data2$long, y=data2$lat, col="slateblue", cex=1, pch=20)
for (j in 1:(length(data1$loc)-1)) { # the number of points on land -1
loc1 <- data1[j,]
loc2 <- data1[j+1,]
inter <- gcIntermediate( # calculate the distance
c(loc1[1,]$long, loc1[1,]$lat), c(loc2[1,]$long, loc2[1,]$lat),
n=500,
addStartEnd=TRUE)
lines(inter, col="tomato", lwd=0.8) # draw the line
}
for (j in 1:(length(data2$loc)-1)) {
loc1 <- data2[j,]
loc2 <- data2[j+1,]
inter <- gcIntermediate(c(loc1[1,]$long, loc1[1,]$lat), c(loc2[1,]$long, loc2[1,]$lat), n=500, addStartEnd=TRUE)
lines(inter, col="slateblue", lwd=0.8)
}
inter <- gcIntermediate(c(dataVI[14,]$long, dataVI[14,]$lat), c(dataVI[15,]$long, dataVI[15,]$lat), n=500, addStartEnd=TRUE)
lines(inter, col="purple", lwd=0.8) # Intersection at Italy and Greece
text(dataVI$loc, x=dataVI$long, y=dataVI$lat, col="black",cex=0.5, pos=3)
Domestic flights network of Delta Airlines
pal <- colorRampPalette(c("#333333", "white", "#00bfff")) #gradient: black-grey-white-blue
colors <- pal(100)
maps::map("world", col="gray17", fill=TRUE, bg="#000000", # dark background
lwd=0.05,xlim=c(-171.738281, -56.601563), ylim=c(12.039321, 71.856229))
fsub <- flights[flights$airline == "DL",]
maxcnt <- max(fsub$cnt) #maximum of the number of flights
fsub <- fsub[order(fsub$cnt),]
for (j in 1:length(fsub$airline)) {
air1 <- airports[airports$iata == fsub[j,]$airport1,]
air2 <- airports[airports$iata == fsub[j,]$airport2,]
inter <- gcIntermediate(c(air1[1,]$long, air1[1,]$lat), c(air2[1,]$long, air2[1,]$lat), n=100, addStartEnd=TRUE)
colindex <- round( (fsub[j,]$cnt / maxcnt) * length(colors) )
lines(inter, col=colors[colindex], lwd=0.8)
}
library(ggmap)
library(curl)
library(httr)
library(dplyr)
Instructions for Google platform is at https://www.littlemissdata.com/blog/maps First, we need to us get_map() function to obtain the information of the address. Then draw the map by ggmap(). The function qmap() is like the combination of get_map and ggmap. It can help us get the map directly.
UClayer_goo_ter <- get_map(location = "University of Chicago", zoom = 15, source = "google", maptype = "terrain")
UClayer_goo_sat <- get_map(location = "University of Chicago", zoom = 15, source = "google", maptype = "hybrid")
UClayer_goo_roa <- get_map(location = "University of Chicago", zoom = 15, source = "google", maptype = "roadmap")
UClayer_sta_ter <- get_map(location = "University of Chicago", zoom = 15, source = "stamen", maptype = "terrain")
UClayer_sta_wat <- get_map(location = "University of Chicago", zoom = 15, source = "stamen", maptype = "watercolor")
UClayer_sta_ton <- get_map(location = "University of Chicago", zoom = 15, source = "stamen", maptype = "toner")
Chicagolayer <- get_googlemap(center = c(-84.50, 39.137580),
zoom = 12, scale = 2,
maptype ='terrain',
color = 'color')
Houstonlayer <- get_map(location = "Houston, TX", zoom = 11, source = "google", maptype = "roadmap")
USA_lon_lat <- c(left = -128, bottom = 23, right = -65, top =52)
USAlayer <- get_stamenmap(USA_lon_lat, zoom = 5)
q1=qmap('Chicago', zoom=12, maptype='roadmap')
library(grid)
library(gridExtra)
g1=ggmap(UClayer_goo_ter,extent = 'device') + ggtitle("Google Maps Terrain")
g2=ggmap(UClayer_goo_sat,extent = 'device') + ggtitle("Google Maps Satellite")
g3=ggmap(UClayer_goo_roa,extent = 'device') + ggtitle("Google Maps Road")
g4=ggmap(UClayer_sta_ter,extent = 'device') + ggtitle("Stamen Maps Terrain")
g5=ggmap(UClayer_sta_wat,extent = 'device') + ggtitle("Stamen Maps Watercolor")
g6=ggmap(UClayer_sta_ton,extent = 'device') + ggtitle("Stamen Maps Toner")
grid.arrange(g1,g2,g3,g4,g5,g6,nrow = 2)
violent_crimes <- subset(crime, offense != "auto theft" & offense != "theft" & offense != "burglary")
violent_crimes$offense <- factor(violent_crimes$offense,
levels = c("robbery",
"aggravated assault",
"rape",
"murder"))
head(violent_crimes)
## time date hour premise offense beat
## 82729 2010-01-01 01:00:00 1/1/2010 0 18A murder 15E30
## 82730 2010-01-01 01:00:00 1/1/2010 0 13R robbery 13D10
## 82731 2010-01-01 01:00:00 1/1/2010 0 20R aggravated assault 16E20
## 82732 2010-01-01 01:00:00 1/1/2010 0 20R aggravated assault 2A30
## 82733 2010-01-01 01:00:00 1/1/2010 0 20A aggravated assault 14D20
## 82757 2010-01-01 02:00:00 1/1/2010 1 13R robbery 2A50
## block street type suffix number month day
## 82729 9600-9699 marlive ln - 1 january friday
## 82730 4700-4799 telephone rd - 1 january friday
## 82731 5000-5099 wickview ln - 1 january friday
## 82732 1000-1099 ashland st - 1 january friday
## 82733 8300-8399 canyon - 1 january friday
## 82757 5200-5299 center st - 1 january friday
## location address lon lat
## 82729 apartment parking lot 9650 marlive ln -95.43739 29.67790
## 82730 road / street / sidewalk 4750 telephone rd -95.29888 29.69171
## 82731 residence / house 5050 wickview ln -95.45586 29.59922
## 82732 residence / house 1050 ashland st -95.40334 29.79024
## 82733 apartment 8350 canyon -95.37791 29.67063
## 82757 road / street / sidewalk 5250 center st -95.41530 29.77119
HoustonCr <- ggmap(Houstonlayer)
HoustonCr +
geom_point(aes(x = lon, y = lat, colour = offense),
data = violent_crimes,
size=1,
alpha=0.2)
HoustonCr +
stat_density2d(aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
size = 2, bins = 4, data = violent_crimes,
geom = "polygon")
HoustonCr +
stat_density2d(aes(x = lon, y = lat, fill = ..level..),
bins = 5, geom = "polygon",data = violent_crimes,alpha=0.4) +
facet_wrap(~ offense)